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ABSTRACT 

The  specific  impetus  for  this  work  was  a  conceptual 
design  of  a  submarine  using  the  toroid  as  the  pressure  hull. 
The  designers  could  not  find  a  ready  body  of  knowledge  to 
obtain  scantlings  for  thier  pressure  hull. 

This  work  began  with  a  review  of  efforts  to  solve 
complete  toroidal  structures.   Several  works  were  found  which 
addressed  general  shells  and  extended  into  partial  toroids, 
but.  the  solution  of  a  complete  toroid  was  not  found  to  be  a 
common  exercise.   Some  of  the  these  works  are  briefly 
reviewed.   An  attempt  was  then  made  to  solve  for  the 
displacements  in  a  thin  walled  circular  toroid  using  the 
energy  method.   Several  problems  were  identified  associated 
with  the  structure  geometry  which  make  the  solution  for  the 
complete  toroid  difficult.   In  addition,  the  functional  used 
for  the  energy  method  needs  to  be  more  complex  than  the 
simple  trigonometric  or  power  series  functionals  used  in  this 
work.   These  two  areas,  geometry  and  functionals,  are  fertile 
areas  for  further  study 
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CHAPTER  1 
INTRODUCTION 

The  purpose  of  this  work  was  to  determine  the  response 
of  a  thin  toroidal  shell  structure.   The  initial  goal  seemed 
achievable  until  the  literature  was  reviewed.   A  simple  and 
understandable  solution  was  not  available.   The  disparity  in 
works  was  apparent  in  such  basic  elements  as  coordinate 
systems.   The  solutions  were  generally  directed  towards 
displacements  or  forces.   The  solution  paths  varied  from 
complex  algebra  to  finite  difference  method.   The  solutions 
were  as  varried  as  the  number  of  individuals  attempting  them. 
It.  was  impossible  to  expect  any  simple  correlation  among  the 
few  people  who  attempted  to  address  this  topic. 

The  result  of  this  work  is  the  identification  of  two 
areas  where  more  time  and  effort  is  required. 

A.  The  impact  of  geometric  curvature  on  the  solution. 

B.  The  impact  of  the  shell  geometry  on  the  assumed 
functional  used  in  the  energy  method  of  the  solution. 

To  be  very  specific,  these  two  areas  have  precluded  this 
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effort  from  achieving  the  initial  purpose  of  this-  year  long 
effort. 

A  word  needs  to  be  addressed  towards  assumptions  that 
are  made.   The  assumptions  are  fairly  standard  for  a  shell 
being  analysed  in  the  linear  elastic  range.   The  assumptions 
are: 

1.  The  shell  is  made  of  isotropic  and  homogenious 
material  which  obeys  Hook's  Law. 

2.  The  thickness  of  the  shell  is  constant. 

3.  The  thickness  of  the  shell  is  small  compared  to 
the  radii  of  curvature. 

4.  A  straight  line  normal  to  the  middle  surface 
before  deformation  remains  straight  and  normal  to  the 
middle  surface  after  deformation  and  retains  its 
original  length. 

In  order  to  establish  a  consistent  presentation, 
geometry  will  be  discussed  here  and,  unless  otherwise  noted, 
all  representations  of  geometry  will  follow  this  coordinate 
system  and  the  definitions.   Solutions  by  other  individuals 
will  be  translated  into  this  reference  system. 
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DEFINITIONS 

a   :     Ratio  of  R  to  r;   a  =  £ 

r 

Since  R  >  r  to  form  a  complete  toroid.  a  will  always 

be  greater  than  1. 
C   :     Circumference  of  the  parallel  circle.:   C  =  2  n   r* 
r   :     Radius  of  curvature  in  the  0  direction  of  the 

unloaded  shell  envelope,  or  radius  of  the  circle 

which  is  rotated  about  the  Z  axis  Cat  distance  R)  to 

form  the  toroid. 
r   :     Radius  to  any  point  on  the  unloaded  toroid  measured 

perpendicularly  to  the  Z  axis. 

ill 

r  =  R  +  r  sin©  =  r  (a  +  sin0). 

R   :     Radius  of  rotation  about  the  Z  axis  of  the  circle  of 
radius  r  used  to  form  the  toroid. 

t   :     Thickness  of  the  plate 

cp        :  Angle  of  rotation  about  the  Z   axis  measured  counter 

clockwise  from  the  positive  X  axis  in  the  X-Y  plane. 

0   :    Angle  of  rotation  measured  clockwise  from  the  local 
perpendicular  to  the  X-Y  plane  in  the  positive  Z 
directon  at  distance  R  from  the  origin  in  any 
direction  <p.      See  Figure  1.   It  sounds  like  a 
lot.  but  it  is  not. 


Figure  i 
Layout  of  Variables 
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p   :    Radius  of  curvature  in  the  <p   direction. 

p  =  — —  +  r 
sin© 

Coordinate  System:  The  standard  right-handed  XYZ  coordinate 
system  is  the  global  reference  system.   Other 
coordinates  (eg.  6  and  <£)  build  on  this  basic 
reference  system. 

Geometric  Curvature:  Curvature  of  the  shell  due  solely  to  the 
geometry  or  curvature  of  the  unloaded  toroid. 

Load  Curvature:  Curvature  of  the  shell  caused  by 

displacements  in  response  to  the  applied  load.  This 
could  also  be  looked  upon  as  the  change  in  curvature 
from  the  reference  (unloaded)  curved  surface. 

Axis  of  Symmetry:  Axis  about  which  the  small  circle  is 
rotated  to  form  the  toroid  -  the  Z   axis. 

Plane  of  Symmetry:  Plane  which  bisects  the  toroid  -  the  X-Y 
plane 

Meridian:  Intersection  of  any  plane  containing  the  Z   axis 

with  the  toroid.   On  either  side  of  the  Z   axis,  the 
intersection  results  in  a  circle  (for  the  unloaded 
toroid)  with  0  as  the  variable.  (See  Figure  2) 
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Parallel  Circle:  Intersection  of  any  plane  parallel  to  the 
X-Y  plane  with  the  toroid.   (See  Figure  2) 
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Figure  2 


Parallel  Circle  and  Meridian 
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CHAPTER  2 
HISTORICAL  VIEW  AWD  THE  PROBLEM 

Historical  View: 

Upon  starting  to  research  the  analysis  of  circular 
toroidal  shells,  the  field  was  found  to  be  not  very  large. 
Books  on  shell  theory  and  analysis  would  neglect  or  only 
briefly  mention  the  complete  toroid.   This  could  imply  that 
the  toroid  was: 

CI)  a  simple  extension   of  general  shell  theory  or 

(2)  much  more  complex  than  general  shell  theory. 
The  latter  seems  to  be  the  case. 

Senjanovic  tl]  gives  a  brief  history  of  shell  theory 
in  his  book.   His  Russian  pride  shows  in  some  instances 
("...shortcomings  were  eleminated  by  representatives  of  the 
Russian  School..."),  but  the  relative  youthfulness  of  the 
field  is  obvious.   Initial  work  in  what  is  now  called  shell 
theory  was  by  H.  Aron  (1874)  and  A.  Love  (1888)  a  mere 
century  ago.   Senjanovic  does  cite  Flugge  -  another  work 
referenced  here  -  as  one  of  the  developers  of  the  classical 
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(or  western)  theory  of  shells.  Yet  even  Flugge  devotes  only  four 
pages  12]  to  the  toroid.  The  following  works  are  illustrative  of 
the  background  on  toroidal  shells. 

V.  Flugge:  Flugge  starts  with  a  generalized  shell  and 
establishes  force  equlibrium.   From  this  equlibrium 
he  solves  differential  equations  to  obtain  solutions 
for  forces  and  displacements.   His  discussion  of  the 
toroid  points  out  some  of  the  problems  inherent  in 
the  full  toroidal  solution. 

L.  Sobel:  Sobel  [3]  was  an  understudy  of  Flugge. 
His  doctorial  dissertation  starts  with  the  same 
generalized  shell  theory  but  is  slanted  towards  the 
buckling  analysis  of  the  toroid.   His  solution  was 
worked  on  a  computer  and  is  difficult  to  follow. 
I.  Senjanovic:  Senjanovic  starts  with  a  force 


ballance  on  a  differential  element  in  a  localy 

orthogonal  coordinate  system,  then  uses  several 

methods  to  solve  them.   Most  notably,  he  employed 

complex  variables  and  some  computer  solutions  circa 

1970. 

S.  Timoshenko:  Timoshenko  [4  3  starts  with  a  plate 
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and  develops  into  a  generalized  shell.   His  attention 
to  the  toroid  is  of  the  same  order  as  Flugge.   This 
writer  found  Timoshenko's  general  development  of 
theory  easier  to  follow  than  others  and  hence  the 
references  to  Timoshenko  in  this  area  may  be  more 
frequent  than  to  other  authors  listed  here. 
E.  Tsui:  Tsui  [53  starts  with  the  constitutive 


relationships  and  develops  them  into  forces.   However 
for  the  toridal  shape  his  solutions  are  on  globally 
orientated  forces  and  displacements.   Tsui  intended 
his  work  more  as  a  handbook  than  a  textbook.   His 
solutions  are  tables  of  values  from  computer  analysis 
of  toroidal  shells.   He  does  provide  however,  a 
consistant  development,  though  shorter  and  in  less 
detail  than  some  others,  for  various  shells  of 
revolution. 
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The  Problem: 

The  specific  problem  to  be  addressed  in  this  thesis  is 
the  response  of  a  toroidal  shell  to  hydrostatic  presure. 
More  specifics  on  the  loading  will  be  addressed  in  chapter 
three.   This  is  to  set  the  stage  for  the  problem  itself 
and  lay  out  the  difficult  areas. 

The  goal  of  any  structural  solution  is  to  take  into 
account  the  loading,  geometry  of  the  structure,  and  the 
material  in  the  structure  and  determine  the  response  of  the 
entire  structure.   The  response  is  then  compared  to  a  desired 
standard  or  a  set  of  failure  criteria  and  adequacy  of  the 
design  is  determined.   This  approach  is  generally  slanted 
towards  displacements  of  the  structure  since  only  through 
displacements  are  the  applied  loads  distributed  and 
equil ibrated. 

In  any  structure,  the  various  elements  are  in 
communication  with  one  another.   The  beam  under  a  simple 
tensile  axial  load  is  the  easiest  example,  but  a  plate  or 
shell  must  also  satisfy  this  connectivity.   Unfortunately, 
the  complexity  of  this  connectivity  or  communication  among 
elements  increases  rapidly.   In  the  axially  loaded  Ctensile) 
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beam,  only  the  axial  displacement  (say  the  X  direction)  is 
important.   Load  the  same  beam  laterally  and  now  the  X  and  Z 
directions  are  required  to  describe  the  beams  response. 
Going  to  a  plate  now  requires  the  X,  Y,  and  Z  directions. 
The  shell  now  takes  the  plate  out  of  the  convenient  XY  plane 
for  a  reference. 

With  appropriate  coordinate  system  selection,  some  very 
practical  shells  can  be  handled  easily.   The  cylinder  and  the 
sphere  are  classical  examples.   The  transition  to  a  toroid 
would  -  by  simple  implication  -  be  very  easy.   Start  with  a 
cylinder  of  radius  r  and  length  of  2nR  and  bend  it  around 
upon  itself.   Discounting  a  few  wrinkles  on  the  inner 
sections  of  the  cylinder,  we  now  have  a  toroid.   The  wrinkles 
however  carry  a  much  deeper  implication. 

The  structure's  ability  to  distribute  the  applied  load 
changes  drastically.   The  complete  toroid  has  no  boundary 
conditions  comparable  to  the  cylinder's  boundary  conditions 
at  X=0  and  X=L.   These  are  replaced  by  connectivity  <or 
continuity)  requirements.   The  analysis  of  toroidal  segments 
is  facilited  by  the  ability  to  insert  boundary  conditions. 
But  even  this  facility  can  be  of  limited  use  if  too  much  of 
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the  structure  is  included  in  the  analysis.  The  simple 
thin  walled  analysis  on  a  cylinder  for  a  force  balance 
becomes  an  algebraic  exercise  of  no  small  feat  for  the 
toroid.   CThis  will  be  shown  in  chapter  five,) 

It  should  be  noted  that  all  of  the  individuals  above 
stated  or  implied  some  difficulty  in  obtaining  a  solution  for 
a  complete  toroid.   The  most  blatant  was  Tsui  [6  J  when  he 
states  that  in  problems  involving  the  complete  toroid,  the 
crowns  C  e  =  0°and  180°  -  see  FIG  1)  can  not  be  solved 
due  to  singularities.   Flugge  [73  gives  a  solution  for  the 
forces  but  states  that  they  cannot  be  used  "...  in  the 
vicinity  of  the  top  and  bottom  circles"  CTsui's  crowns) 
"without  additional  bending. . . ".   He  then  provides  a 
displacement  expression  containing  terms  which  blow  up  at 
0  ■  0°  and  180°. 

With  all  of  this  going  against  an  easy  solution,  two 

courses  of  action  were  undertaken: 

A.  The  singularities  at  the  crowns  could  be  avoided 

by  using  some  expression  which  avoids  the  offensive 

term  C — — )  Csee  Chapter  6).   This  was  avoided  by 
sinO 

referencing  the  differential  element  to: 
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r   =  R  +  r  sin  6  -   r  C  a  +  sin  e  )* 
The  biblical  addage  still  holds:  If  thine  eye  offends 
thee,  pluck  it  out. 

B.  The  method  of  solution  selected  would  be  one  to 
avoid  further  problems  as  pointed  out  by  Flugge  [81. 
The  energy  method  would  avoid  messy  mathematical 
stumbling  blocks  in  the  solution,  since  the  form  of 
the  solution  is  assumed  going  into  the  energy  method. 
This  was  almost  true  as  shall  be  pointed  out  latter. 
Additionally,  this  writer  has  not  seen  this  attempted 
by  others.   Therefore,  when  the  line  forms  to  say 
that  this  method  will  not  work,  it  should  be  a  very 
short  line  indeed. 
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CHAPTER  3 
CONSTITUTIVE  RELATIONSHIPS 

An  expression,  or  series  of  expressions,  relating  the 
material  motion  to  properties  such  as  stress  and  strain  Cand 
further  to  forces  and  moments)  is  required.   This  follows  the 
definition  of  constitutive  equation  of  the  material  as  given 
by  Fung  [9],  "...mathematical  expression  of  the  mechanical 
property  of  a  material. " 

First  we  must  define  our  displacements  (see  Figure  3). 
Three  displacements  will  be  used: 

u  : Displacement  in  the  positive  ©  direction 
v  : Displacement  in  the  positive  0  direction 
w  : Displacement  in  the  direction  of  the  positive 
(ie.  outward  pointing)  normal  at  the  point  of 
interest. 
Rotations  will  be  defined  as  follows  Csee  Figure  3): 

Dq   :  Rotation  of  a  unit  normal  about  the  tangent  to 
the  meridian  in  the  positive  u  direction:  that  is. 
rotation  about  the  positive  ©  axis. 
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Figure  3 


Dispalcements  and  Rotations 
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o  :  Rotation  of  a  unit  normal  about  the  tangent  to 

the  parallel  circle  in  the  positive  v  direction; 

that  is,  about  the  local  positive  <p   axis. 

a     :  Rotation  about  the  positive  w  direction. 
The  loading  will  impact  greatly  on  the  final  outcome. 
For  this  problem,  the  loading  will  be  assumed  to  be  a 
constant  hydrostatic  pressure.   This  has  several 
implications: 

A.  Hydrostatic  loading  always  presents  a  loading 
normal  to  the  surface. 

B.  Because  of  the  global  nature  of  the  hydrostatic 
loading,  any  benifit  to  be  garnered  from  symmetric 
conditions  would  be  applicable  to  this  problem. 

The  second  aspect  of  the  loading,  axisymmetric  loading, 
results  in  the  number  of  unknowns  being  reduced  dramatically. 
Both  Timoshenko  [10]  and  Flugge  [11]  state  that  the  total 
deformation  can  be  represented  by  only  two  displacements;  one 
normal  to  the  surface  and  one  along  the  meridian  in  the 
positive  ©  direction.   Flugge  also  states  [12]  that  that 
stresses  are  independant  of  <p   and  in  fact  any  derivatives 
with  respect  to  <p   are  zero.   Timoshenko  L 1 3 3  agrees,  but 
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the  derivative  going  to  zero  is  implied,  not  stated  as 

bluntly  as  in  Flugge. 

To  take  this  one  step  further,  a  lot  of  things  will  go 

to  zero  because  of  ax i symmetric  loading.   Sheer  stress  and 

strain  are  cross  coordinate  terms,  ie. 

[*]    =  L  (  L  [*]  ) 
e<p       66       6<p 

But  anv  ^  =   0.  so  L*]^  =  l*]^  =  0. 

Because  of  the  type  of  loading  described,  v  and  o  are  not 

required  for  the  description  of  motion.   v  is  not  required 

because  of  the  independance  of  the  motion  on  <f>,    and  o  is  not 

required  because  it  depends  on  a  partial  of  4>.      All  of  this 

is  the  result  of  axisymmetric  loading. 

STRAINS 

Midplane  strains  are  as  follows: 

Cq   :  Is  the  strain  in  the  positive  6   direction,  Sq 
is  made  up  of  two  parts  (see  Figure  4): 

1 .  Cq   due  to  u 

At  two  points  on  the  meridian:  P  and  Q 

Up  =  u     uQ  =  u  +  g§  de 

f  =  £i  =  J_cu  +  £ide-u)  =  L£i 

1   r  dO      60  r  66 

2 .  Cq   due  to  w 
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Figure  4 
Midplane  Strains 
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Length  of  a  segment  is  equal  to  the, radius 

times  the  angular  rotation.   The  original 

length  is  r  d0.   At  point  P  the  radius 

increases  by  w  to  (r  +  w)  so  the  new  length 

would  be  Cr  +  w)  d0.   Similarly  at  point  Q  the 

radius  increases  but  is  now 

(r  +  w  +  —  d0).   The  new  length  would  be 
60 

Cr  +  w  +  £2£d0)d0.   Averaging  the  two  lengths 
60 

at  P  and  Q  gives  an  average  Al  =  wd0  +  j  &=fle- 

Dividing  by  the  original  1  =  r  d0  gives: 

e  -     %     tw  +  -  —  d0]    Neglecting  second 
r       2  60 

order  terms  yields:     c  =  ~r 
Total  Cq   is  therefore: 

y   r  60 
c  -  :  Is  the  strain  in  the  positive  <p   direction 

(see  Figure  4).   Since  the  body  is 

symmetrically  loaded,  all  0's  have  the  same 

displacements.   Therefore  c,   can  be  calculated 

by  calculating  the  change  in  the  circumference 

of  the  parallel  circle  at  <p   =  constant.   Since 

the  circumference  is  related  to  the  radius, 
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the  radius  (r  )  will  be  used  to  determine 


strain  c±. 

* 


C  =  2  n   r 


At  point  P.  displacements  w  and  u  contribute 
to  the  Ar  as  follows: 

|b 

Ar  =  u  cos©  +  w  sin© 

AC  =  2  tt  Ar* 

c,   -   ££  =  L,  [u  cos6  +  w  sine] 
0    C   r 


ROTATIONS 


cj^  is  determined  from  two  points,  P  and  Q,  on  a  meridian 
seperated  by  displacement  u  in  the  positive  e  direction  (see 
Figure  5).   The  unit  normal  vector  at  P  must  rotate 
through  dG  to  go  from  P  to  Q  due  to  displacement  u,  ie. 
r  de  =  u 

.  • .  de  =  i  =  o4 

If  w  was  constant  from  P  to  Q,  no  additional  rotations  would 

be  required  since  the  surface  displaces  uniformly  through  de. 

If  n  ^  0,  then  a  rotation  of  the  local  normal  at  0  will  be 

required  to  maintain  its  normal  relationship.   This  rotation 

wilJ  be: 

d  =  -  SEde/Cr  de>=  -  iSl 
2     60  t66 
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Figure  5 


Rotation  o 
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This  makes  the  total  o  • 

o  =  I   Cu  -  & 
*   r      6G 

Because  of  the  symmetric  loading  and  no  subsequent 

variations  in  the  displacements  with  0,  no  v  displacements  or 

—  exists.   However  a  G&   can  exist  Csee  Figure  6), 
69  ^ 

Due  to  symmetry  in  loading,  |w_|  is  constant  in  <p. 
However,  for  o^  to  change  direction,  an  additional  rotation 
perpendicular  to  o^  must  exist.   This  is  the  rotation  o^. 

Qq=  -  Q^  cos0  d<j> 
Figure  6  gives  a  graphical  presentation  of  Oq. 
Timoshenko  L14  3  gives  a  similar  result  but  achieves  this  as 
a  side  result  of  obtaining  curvature. 
LOAD  CURVATUKES 

Load  curvatures  are  as  follows: 

Xq.    Change  in  curvature  along  the  positive  ©  axis 
due  to  the  applied  load. 

xr.    Change  in  curvature  along  the  positive  <p   axis 
due  to  the  applied  load. 
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Figure  6 


Rotation  o^ 
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It  is  significant,  to  this  writer  at  least,  to  seperate 
these  curvatures  from  the  geometric  curvature  defined  in 
Chapter  1.   The  geometric  curvature  gives  one  an  indication 
of  how  difficult  the  problem  will  be  [15]  and  defines  the 
initial  or  unloaded  reference  surface  for  further  work.   The 
load  curvature  is  the  sole  result  of  the  structure's  response 
to  the  load  applied  and  directly  related  to  the  displacements 
use  to  equilibrate  the  load. 

Load  curvature  along  the  0  direction  is  the  change  in 
rotation  about  the  6  axis  ( — *0  divided  by  the  original 


66 


length  CrdO) 


ka  =  1 —  £-  Co.)  d9 
6   rd9  66       * 

Recalling:  a^   =  ^  (u  -  — ) 
^   r      66 

„  _  1  (6u   _  <52w, 

Timoshenko  [161  gives  a  description  of  this  and  x±   is 
similarly: 

Xx  =  oe  dO/(r*  d8) 

_    cos  6  (u  -  ^w) 

r2Ca  +  sin9)      W 

Ye   can  now  describe  al]  the  pertinent  engineering 
parameters  in  terms  of  displacements.   Adding  to  this  list 
the  general  expression  for  planer  stress,  we  get  the 
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following: 

°e  =  J^a  Cee  +  PV 

=  -^-^  (L  t^S  +  w]  +  £ Cu  cose  +  w  sin63) 

1-a,a   r   <x>        rCa  +  slne> 

o,    =   -£-—   C ^ tu  cose  +  w  sine]  +  £  [^  +  w] 

*"    1-v2   rCa  +  sine)  r   <5e 

This  almost  completes  our  constitutive  relationships. 
Stress,  strain,  rotation,  and  local  curvature  can  be 
expressed  in  terms  of  the  displacements  u  and  *r.   The  last 
step  is  to  proceed  to  normal  forces  and  moments.   This  is  as 

follows: 

h 

;  v  «*  e,  0 


h 

N 

=  /  2 

-h 

2 

h 

O 

dz 

M 

V 

=     /     2 

-h 

O 

z   dz 
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CHAPTER  4 
EQUILIBRIUM  OF  FORCES 


The  next  aspect  of  the  solution  to  be  presented  will  be 
that  of  establishing  equilibrium  conditions  for  the 
differential  element.   The  differential  element  is  shown  in 
Figure  7.   Geometry  rears  it's  ugly  head  in  that  the 
element  needs  two  reference  points  to  describe  it's  surface. 
The  first  is  the  axis  of  rotation  of  the  shell  in  the  & 
direction.   The  second  is  the  origin  of  the  global  coordinate 
system.   The  first  axis  can  be  referenced  to  the  global 
coordinate  system  by  R  and  4>. 

Solutions  for  partial  toroidal  shapes  all  take  advantage 
of  the  known  boundary  conditions  on  the  edges.   When  a 
complete  toroid  is  analyzed,  the  boundary  conditions  double 
back  upon  themselves  and  become  internal  vice  external 
forces.   In  these  instances,  symmetric  loading  and 
deformations  are  invoked.   Tsui  117]  gives  tabulated 
computer  results  for  both  solutions  and  the  differences  are 
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Figure  7 


Differential  Elenent 
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obvious. 

In  the  following  discussion,  and  the  associated  figures, 

the  following  conventions  will  be  used  (see  Figure  6>: 

> 
R  is  a  generalized  force  ie.  a  force  or  moment  in 

the  v  direction  d  «  6,  <t>,    n>.   The  units  for  F  are 

[Force/unit  length]. 

F   is  the  force,  ie. 

F.  =  [Force/unit  length]*[ length]  ■  CForce] 

N6'  N<2>  Normal  forces  in  the  <;:>  direction  acting  on 

the  side  of  the  element  on  which  <£>  is  constant. 

Q^  Sheer  force  acting  on  the  $   ■  constant  face  of 

the  element.   Note:  Due  to  symmetry  of  loading, 

Ma,  M,   Moment  to  induce  rotation  in  the  <]l> 
direction.   The  moment  is  defined  as  positive  with 
respect  to  the  right  hand  rule  as  applied  about  the 
local  positive  <2>  direction. 
Each  force  will  be  addressed  graphically  with  three 
isometric  diagrams  and  a  geometric  breakdown  to  the 
appropriate  F  .   The  F   for  each  force  (N^,  N,,  .  .  ect.)  can 
then  be  summed  for  a  force  balance. 
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Force  Figure 
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Figure  8 


Forces  Defined 
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Figure  P 
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6N, 


F^:  -Nrt  sinp  +  CNA  +  — ^d6>  cosa  sinp 


0'   "6 


'6 


66 
6N 


Fe:  -Ne  cos?  +  CNe  + 
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Figure  10 
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Figure  11 
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F^:   Q@  cosa  sin^  +  Qa  cosa  sin/3 


Contribution  to  F^  doe  to  Q^: 

>  i 

F^:    Qq  sina  cos/3  -   Q^  cosa     sin/3 
F    :    Qq  cosa   cos/3  +   Q^  cosa      cos/3 
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Figure   12 
P 


0   - 


a  ■    \    C    d6    ♦    o,    + 

2  <p 


So 


6) 


<5o, 


2         B    <5«!>  r) 

NOTE:  Z  =  P  Area  =  P  [<R  +  r  sin6)d0]  [r  d6] 


9 


Contribution   to  F      doe  to  Z: 


F.:    -   Z   cosa  sin/3 


F~:     -   Z   cos/3  sine 


F    :    -   Z   cos/3  cosa 


-   39    - 


Figure   13 
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Figure    14 
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Relationships  between  the  N  and  N   must  now  be 


established. 

Ne  e  Ne  r*  de 

NA  +  — 2  =  CN  A+  — H)Cr*  +  2L_>  d<A 

■  NAr  d0  +  NAC^-)d«^»  +  C — £>r  d<*>  +  < — 2)<2£_)d^ 
w  6G  66  66    66 

6Ne    „,  s  m 

~   (NA  +  — — )  r   d0,   neglecting  terms  with— —  as 
66  6e 

being  relatively  small 

Similarly: 

N ,    =  N *  r  d6 

Qe  -  Q'e  r*  d<> 

6ge      ,   6g' 

U     60  ^         SO 

ne   =  Me  r  d* 

6MA       ,    <5M«    ^ 
M.  =  m'  r  d6 
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CHAPTER  5 
ENERGY  METHOD 

Any  attempted  solution  must  maintain  some  sort  of 
equlibrium.   In  the  energy  method,  the  equlibrium  comes  from 
minimizing  the  energy  from  the  external  load  and  the  internal 
energy  of  the  shell.   For  the  external  energy,  a  surface 
integral  of  the  hydrostatic  load  and  the  normal  displacements 
gives  this  value.   The  internal  energy  is  determined  by 
integrating  stress  times  strain  throughout  the  shell  volume. 
The  external  energy  is  subtracted  from  the  internal  energy 
and  the  resulting  functional  is  generally  refered  to  as  n, 
the  total  potential  energy  of  the  structure.   Vith  V 
representing  the  internal  energy  and  W  representing  the 
external  energy,  the  functional  looks  like: 

n  *  V  -  V 

If  the  value  of  n  is  at  a  minimum  value,  the  structure 
wiil  be  in  equlibrium.   This  is  the  principle  of  total 
potential  energy  [183.   This  work  will  not  address  buckling 
or  post  buckling  behavior,  only  the  linear  elastic  region. 
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For  this  situation,  equlibrium  is  achieved  when: 

6n      _  6V_  _  6V_  m   Q 
dxi  ix\  6q 

X  \.  1. 

where  q   is  any  generalized  parameter  used  to  determine  V  and 
V. 

Unfortunately,  the  parameters  used  to  determine  V  and  V 
are  u  and  w  -  which  is  what  we  are  attempting  to  solve.   To 
get  around  this  difficulity,  a  solution  will  be  assumed. 
Then  an  energy  balance  will  be  conducted. 

The  assumption  of  a  solution  is  not  without  risk. 
First,  the  solution  must  meet  known  natural  boundary 
conditions.   Then  the  accuracy  of  the  solution  is  limited  by 
the  closeness  of  the  assumed  solution  to  the  actual  Cand 
unknown)  solution. 

In  the  case  of  the  toroid,  there  are  no  "boundary 
conditions"  as  was  pointed  out  in  chapter  four.   In  place  of 
the  boundary  conditions,  we  have  conditions  of  compatability 
to  limit  the  solution.   As  stated  earlier  CtlP]  and  [20]), 
symmetry  of  deformation  can  be  assumed  for  the  hydrostatic 
loading.   The  plane  of  symmetry  is  the  XY  plane. 


-  44  - 


To  determine  the  internal  energy  we  proceed  as  follows: 

V  ■  I       O  C  dV   ;    i  6  j  ■  6,  4> 

-  J  /^  CN    c    +  M    *   >  dS 

Recall  that,  due  to  the  type  of  loading, 


cross  terms  <v  *   j )  equal  zero. 

v  -  i  ss   CNe  ce  +  %  e0  +  Me  **  +  M*  V  ds 

Eq  1:    N0  ee   =  K  C*6  +  P  c0>  ^     ;   K  =  -^ 
=  K  <£§  +  2  v  e0  e^   +  e£> 

Simi larly: 

EfJ  2:       n0  c0  -  K  ^0  *  2  *  <>  ce  +  eP 

Eq  3:     Me  *e  =  D  C*e  +  *f    *&  ;       D  =  J<^y2>  "  *  <nT> 

=  K  cljjj)  <«g  +  2  P  *e  *^  +  *£> 

^  4:    M0  %  =  K  <rZ>  <*0  +  2  V   *<*>   *e  +  *6> 
Recall  from  chapter  four: 

Eq  5,6:   rA  =  i  C—  +  w) ;  c,  ■  - Cu  cosO  +  w  sine) 

*        r   66         *    rCa  +  sine) 

Eq  7,8:    «A  =  L-  C*H  -  ^) ;     «,  =   -   COs6 Cu  -  <*) 

r2   60    662       *    r2Ca  +  sine)       SO 

The  unknowns  u  and  w  still  persist.   The  assumption  will 
be  made  that  u  and  w  are  functions  of  e.   A  series  solution 
will  be  assumed  as  follows: 
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w(6)  =  I  w  sinCmO),   limits  of  summation  will  be 

addressed  later 
u<6)  =  I  u   —  CsinCmO)) 
■  Z  u   m  cosCmO) 

m 

This  turns  out  to  be  not  a  very  good  assumption  for  the 
the  series  representation  of  the  displacements.   There  are 
some  inherent  limitations;  for  instance  sinCm©>  is  always 
zero  at  0=0.   While  not  in  conflict  with  the  conditions  of 
compatabi lity,  this  is  an  additional  constraint  on  the 
solution  displacements.   If  the  actual  displacements  are  not 
zero  at  0   =  0,  then  there  is  a  guaranteed  error  in  the 
assumed  solution.   The  question  to  be  answered,  as  with  all 
assumed  solutions,  is  how  good  is  the  answer  obtained 
compared  to  the  answer  that  is  required.   More  on  this 
subject  in  Chapter  6. 

For  ease  of  notation,  let: 
<t>     ■  sinCmO) 

m 
> 

0  =  m  cosCm6) 

n-i 


<t>     -   -m2  sinCmO) 
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Therefore: 

Eq  9:         w<0>  ■  Z  w  d> 
Eq  10:        uce>  =  Z  u   <£ 

m   m 

Boundary  conditions  for  this  problem  are  actually 

compatibility  conditions.   At  6  ■  90°  and  270°  the  uC6> 

displacement  should  be  zero  and  — Cw<6>)  should  be  zero. 

SO 

^s^Lciw   0)=lw   — C0  )  -  E  w   ^' 

jf|Q      c'^n       mm  m   **/\    m  mm 

Since  w   are  merely  coefficients  now  in  the  series  solution 

rr. 

they  are  unaffected  by  the  derivative.   To  look  at  this  in 

more  detail: 

—  =  w   cos(O)  +  w   2  cosC26)  +  *#   3  cos<36)  ♦  .  .  . 
66     *  2  3 

+  w   n  cosCnO) 

r, 

For  coiiveince.  let  n  -  10  for  the  trial  solution.   At  6  ■  90° 
and  270°.  cosCmO>  ■  0  for  m  *  1,  3,  5,  .  .   and  cosCm6>  ■  ±  1 
for  m  =  2,  4,  6.  .  .  .   For  the  series  to  be  equal  to  zero 
set  w   =0  for  m  =  2,  4,  6,  .  .  .   Similar  reasoning  gives 
u   =0  for  m  *  2,  4,  6,  .  .  .   CThis  explains  why  the 
functional  for  u   was  selected  as  the  derivative  of  the 

m 

functional  for  w  . )   As  shall  be  seen  later,  the  method  is 

m 

smart  enough  to  realize  this  limitation.) 

The  elemental  area  CdS)  must  be  investigated  for  the 
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toroid.   dS  =  dle  dl  .  :   dle  ■  r  d6   ;  dl  ,  -  <R  4'   r  sin6>  dtf» 

B  r  de  (R  +  r  sine)  d# 

As  a  check:   S  -  f*nf*nr   <R  +  r  sine)  d0  d* 

o   o  ^ 

«  4  n2  R  r.   This  agrees  with  the  CRC 
Standard  Mathematical  Tables  [21]  for  the  surface  area  of  a 

toroid,   therefore  dS  is  correct. 

jcrt 

The  ultimate  objective  is  to  get  21i_  b  o.   In  this  case 

q   =  u   and  w  .   So  for  each  specific  u   -  that  is  u,  -  one 
gets  an  expressison  such  that: 


6TT_  m   6V_  _   SV_  u 
6u ,    6u ,    <5u , 

k        k        k 

Likewise  a  series  of  m  expressions: 

SU      _  6V_  _  6V_  m 
6w,    6w,    <5w, 

k        k        k 


0 


0 


For  a  nominal  m  =  10.  this  yields  20  equations.   Fortunately, 
u   and  w   are  20  unknowns. 

rr,  m 

Now  the  internal  energy  can  be  formated.   Equations  9 
and  10  are  substituted  into  equations  5  through  8.   These 
equations  are  then  inserted  into  equations  1  through  4.   This 
is  combined  with  the  expression  for  dS  to  give  V,  the 
internal  energy.   To  facilitate  the  future  derivatives  which 
will  be  required,  the  subscripts  m  and  n  will  be  used  to 
indicate  the  seperate  terms  in  a  product,  ie. : 
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£%   =  £A  £&   -    i:   CZ  u  4> "  +  Z  w   ^  )  i  (I  u   ^ "  +  Z  W   #  ) 
All  or  the  above  yields: 
V  =  n   r2  K  /**  rkZu  <>"  +  Zw  0  )^CZu  <f>"   +  Zw  <*>  ) 

O    *-r      mm        m  m  r      nn        n  n 


1  » 

x  tW  a  -i  c,m^(cose  Zu  0  +  sin©  Z  w  <f>   ) 

+  C&-)  -4  <Zu  0"  +  Zw  </>")  -4-CZu  <t>"   +  Zw  0"> 

1Z     2       mTm        mrm     2     n  n        nTn 

+  cJ^O  ^  CZu  <p'    +  Zw  0")  — ^^^ CZu  0'  +  Zw  «/> 

1Z   r2     m   m      m  ^  r2Ca  +  sine)    n  n      n  n 


h\      cos6 

r2(a  +  sine) 


c'J-t)  -t; — ^^ CZu  0  +  Zw  4>   > 

1^    „2/_   _.   _i_^-w>»     mm 


r2(a  +  sine)    n  n      n  n 


Next,  carrying  out  the  multiplications  yields: 
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V   =   n   r2    K   S2n\j—  ZZCu   u      0*V   +   u   w      0m0n   *  w   u      0  0** 

«J      ••     z                    m     n         m     m                mnmn  mnmn 

+    W    W        0    0) 

m     n         m     n 

+— — zz<u  u    00'  cose  +  u  w    00    sine 

2,  ^             ^j  ^/>,                m     n         m     m  mnmn 


r    Ca   +   sine) 

+   w  u      00     cose  +  w  w     00     sine) 

iz<u  u     0*0'   cos2e 


r2(a   +   sine)2 


m     n      '  m '  n  m     n      '  m '  n 


+  u  w  00  cose  sine  +  w  u  00  cose  sine 


+   w  w     00     sin2e) 

m     n         m     n 


.2  ..      ..  ..      ..  ..      .. 

+     <■  )Z2.<U     U        00       -UW       00       -WU        00 

A  ^,  „  4-  m     n         m     m  mnmn  mnmn 

12r 

+     W    W        0    0) 

m     n         m     n 

+    c-l^-)      2    *'  cos9  ZSXu   u      0*V    -   u   w      00" 

12r4         Ca     +     Sine)  mnmm  mnmn 

•  ■      >  "      > 

-      W     U  00         +      W     W         00) 

m     n         m     n  mn         mn 

+    (_D— _)    cos   S — _Cu   u00      -uw00 

12r4      Ca   +   sine)2      m    n   m   m  mnmn 

-wu      00     +  w  w     00)lCa  +  sine)    de 

mn      mn  mn         mnJ 

This  is  all  well  and  good  but  actually  - —  is  what  we 

really  need.   Again,  for  this  problem  q   =  u   and  w  j  and 
k  =  1,  2,  .  .  .  10.   Applying  term  by  term  differentiation  is 
straightforward  but  tricky  and  tedious.   A  couple  of 
examples: 

—   Cw   u   0   0)  =  w   0   0' 

j-         m    nmn        mmn 
k 
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But   observe: 

jp  ••  ••  ••  ••  ••  •• 

(um    u      d>      0    )    =    u„   0      0,     +    u      0,0 

=     2     U        0        0, 

m     r  m         k 

Since  m  and  n  can  both  serve  as  a  counter  index,  the  two 
terms  may  be  combined.   The  key  to  this  is  that  both 
functionals  C0  )  have  to  be  of  the  same  order  of  the 
derivative.   Observe: 

Cu   u   0   0  )  =  u   0   0,  +  u   0,  0 

x...       ir<        ri'm'n        m'm'k       nkrr> 
6Uk 

Doing  all  of  the  above  and  collection  like  terms  gives  the 

following  two  expressions: 

The  partial  of  V  with  respect  to  u  . 

6u^  m   °     Tn  k    Ca  +  sine)   k  m  m   k 

2  cos20    ,'   ' 
0„  0i 


Ca  +  sine)2   m   k 
+  ciLr)C2  tf>"*f>"   p»^cose)  (j^'   sine  +  h'  cose) 

12  T m^k      Ca     +     SinWJ     ^k^m  ^rrTk 

+  — 2  cos  e   0*0*)>  ca  +  sine)  de 
ca  +  sine)2   m  k 

+!%»  /2Tr  C2  0'V  +  2_± C*'V  sine 

m   °      k  m   Ca  +  sine)   k  m 

+  ^'  cose)  +  2  cose  sine  0> 
m  k         Ca  +  sine)2   k  m 

+    C^C-2    0,0   -   £  ^  cosu  C0,  0   +  0  0,  ) 
17      k  m   Ca  ♦  sine)   k  m    "■  k 

0'0'))Ca  +  sine)  de] 


Ca  +  sine)2   k  m 
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The   partial    of    V   with   respect   to   w    . 

Q—  -   n   K    pu      /2Tr<2    ^'V     +    2— ^ C^'Vi.    sine  +   <p  ri      cose) 

k 

+  2  cose  sine  .' . 
Ca  +  sine)2   ""  k 

-  cJU)C2  00  ?  ^^uvC^^''  +  4> V !"> 

+  — 2  CQS  e   0 '*S>  Ca  +  sine)  de 
Ca  +  sine)2   m  k 

+Zw   /2"  C2  *0   +   4  "  slne  0  ,0 

+  2  sln2° 0  0 

Ca  +  sine)2   k  m 

+   Cr7)C20  0    +   C^t^m   +   ^m^lr5 

1^      km     Ca  +  sine)    km      m  k 

+  — 2  cos  °  5  0,0>))Ca  +  sine)  de] 
Ca  +  sine)2   k  m  J 

The  external  energy  CV)  should  be  looked  at  next.   The 
energy  is  the  product  of  the  load  -  hydrostatic  pressure 
CP)  in  this  case  which  is  always  normal  to  the  surface  - 
times  the  displacement  along  the  line  of  action  of  the  force 
-  w<e)  -  in  this  case. 

V  ■  /   P  w  ce)  dS 

m 

■  'V  fV*   P£w  0  r  (R+r  sine)  de  d0 

O     O  m  ^m  ^ 

Substitute  R  *  r  a 

-   2   n  P   r2   Z  w     /2Tt  <t>     Ca   +  sine)   de 

m         O         T  m 
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And  the  partials  are: 

$¥—  -   2   n  P  r2   /2Tr  0,     <a  +  sin6>   d6 
6w,  °        k 

k 

**L  =  o 

The  series  expressions  can  now  be  rearranged: 

&YI     m   6V   _  <5W  m   0 
6q.    6qt    6qv 

6V     6W 

- —  =  ^—    >      <?  ■  uu  >   **u 

6q     6q         v     k    k 
This  can  be  represented  in  matrix  format  as  follows: 
Eq  11        ^—  =  [D3  <u>  +  [E]  <w>  b  0  -  6W 


6u ,  6u , 

k  k 

Where  [[)]  represents  the  coefficients  of  u   in  ' 


*■* 


(see  page  9). 

xv 
LE]  represents  the  coefficients  of  w.  in  — — 

k 

(see  page  9) 

<u>T  =  u^  =  Ul,  u2,  u3,      u1Q 

<W>T  ■  Wm  =  *V   *V   *V   '   •   "lO 

Similarly  in  matrix  format: 

Eq  12        ^—  =  [F]  <u>  +  [G]  <w>  =  <P>  -  6V 


6w,  6w, 

k  k 

XV 

Where  [F]  represents  the  coefficients  of  u   in  - — 

Csee  page  10). 

<5V 

[G]  represents  the  coefficients  of  w,  in  — — 

Csee  page  10) 
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<P>T  a  £*_  =  £¥_  +  £*_  +  +6V 

6wk    6*1         6«2         '  6wio 

Solving  equation  11  for  <u>  gives: 


<u>  =  -  [D]"1  CE3  <w> 


Use  this  expression  in  equation  12: 


<[F]  tDl-1  [El  +  [G3>  <w>  *  <P> 


Solve  for  <w>: 


<w>  =  <[F]  [Dl_1  [El  ♦  [Gl>  *  <P> 
If  the  solution  for  w   and  u   are  correct,  then  these 

mm 

can  be  used  to  reconstruct  w<6)  and  uC6)  and  from  these, 
using  the  constitutive  relationships  in  chapter  four, 
stresses  in  the  0  and  <p   directions  can  be  calculated.   Roark 
[221  gives  expected  values  which  can  be  used  for  comparison. 

Some  justification  for  Roark 's  results  should  be 
provided  since  this  is  how  the  solution  will  be  checked.   He 
provides  no  specific  reference  for  these  answers  in  his  book, 
but  working  backwards  one  can  see  that  a  simple  thin  walled 
approximation  analysis  was  performed,  similar  to  that  of  a 
cylinder.   Roark  provides: 

P*     „    ?*   c   a  *   2Sin6  , 

O  .  =  ;     Oa  ■  C  J 

9        2t       "       t     a  +   sine 
(See  Figure  15) 
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Figure  IS 


Thin  Hailed  Analysis 
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o,  is  a  simple  force  balance  in  which  the  pressure  times  the 
internal  cross  sectional  area  equals  the  stress  times  the 
material  area. 

PA  -   a .    A 
P  <n  r2>  ■  o,    C2  w  r>  t 

*        2t 
A  similar  approach  is  used  on  o^.   The  pressure  acts  on  the 

area  A.  which  in  this  case  is  C2  r)C2  n   R) .  Oq   acts  at  6  and 

0  +  180°.   The  thickness  is  constant  at  t  but  the 

circumference  at  O  and  0  +  180°  are  different. 

CCO)  =  2  n   r*  «  2  n   r<a  +  sin©) 
The  metal  area  at  O  is  therefore  AC6)  =  C<0)  t.   The  force 
balance  is: 

p  a  =  oce  )  cce  >  t  +o<e  >  cce^>  t 

11  2        2 

with  6  ■  e  and  e„  =  6  +  180°.   Substitute  Roark's  oC©)  into 

1  2 

the  above  equation  gives: 

Pr   a  +  lsinei 
p  4  n   r  R  =  —  C 2 L)  2  n  rCa  +  sine  )  t 

t   a  +  sin61 

Pr   a  *  lsine> 
+  Li.  c 2 £)  2  n  rCa  +  sine  >t 

t   a  +  sine 

2 

Cancelling  some  terms  shows  the  equality: 

2a  ■  2a    Q.  E.  D. 
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This  must  have  been  how  Roark  arrived  at  these  relationships. 
This  is  how  the  solution  will  be  checked. 

THE  SOLUTION 

The  method  outlined  above  requires  the  integration  of 
several  quantities  which  are  not  readily  found  in  a  table  of 
integrals.   The  integration  was  done  using  Simpsons  Rule 
[233.   It  was  found,  using  known  integrals,  that  50  steps 
with  in  the  interval  0  to  2rt   gave  good  accuracy  C5  to  6 
significant  digits  accuracy).   A  program  was  written  to 
accomplish  the  integration  and  the  subsequent  matrix 
manipulations  to  solve  for  <u>  and  iw>. 

To  check  the  output  from  the  program,  a  spreadsheet  was 
set  up  to  handle: 

%Ke>  =  2  m     4,  6wce)  m  ^  w    ^ 

u<e>  =  z  u     0*  6uce»  m  z  u     0- 

it*    Tin  j-/\  m    m 

The  spread  sheet  also  calculated  Cq,    €*,    ovx,  and  o>.       As   a 
check,  o^  was  averaged  and  used  to  calculate  a  pressure  to 
compare  with  the  input  pressure. 
THE  RESULTS 

The  results  are  disheartening.   In  the  sample  output 
provided,  only  one  of  the  stress  values  used  as  a  comparison 
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came  close  to  the  expected  value.   a,  should  have  been 

constant  and  have  a  value  of  5000  psi.   It  was  not  constant 

but  it  did  manage  to  equlibrate  most  of  the  applied  pressure 

C99.9%).   This  is  pretty  good  for  an  approximation.   The  Oq 

is  another  story. 

Oq   is  too  low  and  does  not  follow  the  solution  in  Roark 

P  r  c   a  *  lsine  , 
t    a  +  sine 


The  following  are  the  default  values  in  the  program: 
GEOMETRY 


Radius  of  rotation 

Radius  of  the  circle 

Thickness  of  the  shell 

MATERIAL 

Poisson's  ratio 

:  3E+07 

SOLUTION 

Number  of  steps  in  the  Simpsons 

integral 
used  in 

the  series  approximation 
Pressure 


8  inches 
2  inches 
. 02   inches 

. 3  Modulous  of  elasticity 


50  Number  of  terms  to  be 

:  10 
:-100   psi 
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The  program  output  for  w  and  u   are: 


m 
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 


-1.213117E-03 

6.043412E-09 
-1 . 034239E-04 

3.250964E-09 
-3.933767E-05 

3. 328626E-09 
-1. 16143E-05 

3. 831875E-09 
-2.034983E-06 

2.81034E-09 


u 

m 

-1.137496E-03 

1.491326E-09 
-1.151281E-05 

2.026468E-10 
-1.579169E-06 

9.227051E-11 
-2.382749E-07 

5.9864E-11 
-2.535156E-08 

2.817199E-11 


Notice  that  the  even  numbered  terms,  which  are  to  be  set 
to  zero  as  indicated  earlier,  are  several  orderes  of 
magnitude  below  most  of  the  odd  values. 

The  reconstructed  output  from  the  program  is  as 

follows: 
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The  spreadsheet  out  put  to  check  the  solution  is: 


P  = 

-100.00 

psi 

-99.90  < 

Sin  Pres 

s  Cpsi) 

V    s 

0.30 

E 

3.3E+07 

E  = 

3E+07 

K   = 

8.00 

in   a  ■ 

4.00 

r  = 

2.00 

in 

t  = 

0.02 

in 

e 

w 

6w 
60 

u 

6u 

6© 

°e 

°<t> 

C°) 

Cin) 

rad 

Cin) 

rad 

Cpsi) 

Cpsi) 

0 

0. 

0E+00 

-1. 

8E-03  ■ 

-1.2E-03 

0. 

0E+00 

-1461 

-4870 

15 

-4. 

4E-04 

-1. 

4E-03 

-1.1E-03 

4. 

2E-04 

-1710 

-4738 

30 

_r7 

2E-04 

-8. 

1E-04 

-9.8E-04 

6. 

8E-04 

-1946 

-4607 

45 

-9. 

OE-04 

-5. 

7E-04  ■ 

-7.8E-04 

8. 

4E-04 

-2123 

-4405 

60 

-1. 

OE-03 

-4. 

2E-04  ■ 

-5.4E-04 

9. 

6E-04 

-2258 

-4248 

75 

-1. 

1E-03 

-2. 

2E-04  ■ 

-2.8E-04 

1. 

OE-03 

-2345 

-4159 

90 

-1  . 

1E-03 

-2. 

3E-08 

-2. 8E-08 

1. 

1E-03 

-2375 

-4131 

105 

-1. 

1E-03 

2. 

2E-04 

2.8E-04 

1. 

OE-03 

-2345 

-4159 

120 

-1. 

OE-03 

4. 

2E-04 

5.4E-04 

9. 

6E-04 

-2258 

-4248 

135 

-9. 

OE-04 

5. 

7E-04 

7.8E-04 

8 

4E-04 

-2123 

-4405 

150 

-7 . 

2E-04 

8. 

1E-04 

9.8E-04 

6, 

8E-04 

-1946 

-4607 

165 

-4. 

4E-04 

1. 

4E-03 

1.1E-03 

4. 

2E-04 

-1710 

-4738 

180 

-9. 

6E-08 

1. 

8E-03 

1 . 2E-03 

9. 

2E-08 

-1461 

-4870 

195 

4 

4E-04 

1. 

4E-03 

1.1E-03 

-4. 

2E-04 

-1269 

-5190 

210 

7. 

2E-04 

8. 

1E-04 

9.8E-04 

-6. 

8E-04 

-1085 

-5498 

225 

9 

OE-04 

5. 

7E-04 

7.8E-04 

-8. 

4E-04 

-894 

-5654 

240 

1 

OE-03 

4. 

2E-04 

5.4E-04 

-9, 

6E-04 

-747 

-5769 

255 

1 

1E-03 

2, 

2E-04 

2.8E-04 

-1 

OE-03 

-659 

-5854 

270 

1 

. 1E-03 

6. 

8E-08 

8.4E-08 

-1 

1E-03 

-630 

-5886 

285 

1 

. 1E-03 

-2. 

2E-04 

-2.8E-04 

-1 

. OE-03 

-659 

-5854 

300 

1 

OE-03 

-4. 

. 2E-04 

-5. 4E-04 

-9, 

. 6E-04 

-747 

-5769 

315 

9 

, OE-04 

-5 

. 7E-04 

-7.8E-04 

-8 

. 4E-04 

-894 

-5654 

330 

7 

. 2E-04 

-8 

. 1E-04 

-9.8E-04 

-6 

. 8E-04 

-1085 

-5498 

345 

4 

. 4E-04 

-1 

. 4E-03 

-1.1E-03 

-4 

2E-04 

-1269 

-5190 

360 

1 

. 9E-07 

-1 

. 8E-03 

-1.2E-03 

-1 

. 8E-07 

-1461 

-4870 

Note:   Due  to  the  precision  of  the  machine  being  used,  the 
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forced  compatability  conditions  appear  to  be  non-zero.   The 
forced  conditions  at  0  =  90°  and  270°  are  significantly  less 
than  the  rest  of  the  solution  and  therefore  satisfactory. 
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CHAPTER  6 
CONCLUSIONS 

The  apparent  simplicity  of  the  toroidal  shell  belies  the 
complexity  of  the  mathematics  to  analyse  this  shell 
structure.   The  analysis  of  simpler  shells  Cplates, 
cylinders,  and  spheres)  has  the  advantage  of  the  geometric 
curvature  remaining  either  constant  or  at  lease  maintaining 
its  sign  positive  or  negative.   The  geometric  curvature  of 
the  toroid  changes  from  positive  to  negative  as  6  goes  from  0 
to  2h. 

The  impact  in  the  change  of  sign  in  the  radius 
of  curvature  about  the  2  axis  has  been  seen  in  the  results 
of  this  work  and  in  applying  other's  solutions.   This  impact 
is  emperical  and  could  be  the  subject  for  further 
investigation.   The  first  blatant  statement  was  found  in  Tsui 
[24  1  where  he  states,  "...the  relevent  differential  equations 
to  assume  singular  solutions,  and  consequently,  problems 
involving  the  crowns  such  as  a  complete  toroid  cannot  be 
solved. "   The  culpret  was  writing  the  differential  equation 
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in  terms  of  (Tsui's  nomenclature): 

r  =  b  (1  +  -^ — ) 
2  sinO 

(This  papers  nomenclature): 

p=  r  (1  +  Cjp/sine) 

To  avoid  the   — - term,  Tsui  shifts  variables  -  both 

sine 

dependant  and  independant  -  and  provides  influence 
coefficient  matrices.   (These  are  outputs  from  numerical 
analysis. ) 

As  mentioned  earlier,  Tsui  addresses  the  significant 
differences  between  the  partial  toroidal  shape  and  the 
complete  toroid.   Others  may  have  addressed  partial  toroidal 
shapes  tangential ly  but  Tsui  once  again  is  very  blunt  in  his 
statements.   Again  remember  Tsui  is  intending  his  work  as  a 
practical  handbook.   His  intended  audiance  wants  practical 
answers  to  real  problems.   Others,  Flugge  [253 
specifically,  attempt  to  give  solutions  but  he  adds,  "...this 
solution  cannot  be  realized. .. because  it  again  leads  to  an 
incompatability  of  deformations. . . ". 
CURVATURE 

In  attempting  to  point  the  finger  at  why  the  solution 
is  so  difficult,  the  answer  comes  back  to  curvature.   This 
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work  has  seperated  curvature  into  geometric  (.unloaded)  and 

loaded  curvature.   This  writer  feels  that  the  geometric 

curvature  is  the  culpret.   Notice  the  expression  for  the 

radius  of  curvature  in  the  0  direction  given  in  chapter  one: 

p  ■  — —  +  r 
sine 

For  the  range  0°  <  e  <  180°;  sine  >  0  and  p  >  0.   Conversely, 
for  180°  <  O   <  360°;  sine  <  0  and  p  <  0.   At  e  -  0°  and  180° 
sine  =  0  and  p   «  ±  oo.   Two  things  were  observed. 

A.  In  previous  solutions  obtained  in  this  study-  ie. 
worse  than  the  one  presented  -  the  stress  outputs  for 
the  two  regions  differed  markedly  between  0°  <  6  < 
180°  and         180°  <  e  <  360°.   This  could  have 
been  due  to  a  lot  of  reasons,  but  this  was  also 
coincident  with  . . . 

B.  Flugge's  solution  for  the  u  displacement  [26] 
was  being  investigated.   This  tended  to  blow  up  from 
e  *  180°  to  360°/0°  due  to  a  term  lnCtanfo.   Other 
problem  terms  existed  Ccote)  which  agrivated  the 
solution  at  the  points  O   ■  0°  and  180°.   In  just 
"playing"  with  the  solution,  the  points  of  curvature 
change  kept  comming  back  as  critical  points  requiring 

-  64  - 


further  attention. 

Something  happened  at  6  «  0°  and  180°.   In  a  cylinder 
these  two  points  are  no  trouble  -  but  the  radii  of  curvature 
are  «>  and  r  -  and  are  constant!   In  a  sphere  these  two  points 
are  no  problem  -  but  the  radius  of  curvature  is  r  in  all 
directions  -  and  constant!   In  the  toroid  the  radii  of 
curvature  at  both  points  are  p  ■  ±  oo  and  r.   The  key 
difference  is  that  p  is  changing  from  p  >  0  to  p  <  0  and  at 
0   =  0°  and  180°,  *.  Cor  =-)  just  happens  to  be  passing  through 
zero. 

It  was  mentioned  earlier  in  Chapter  4  that  it  required 
two  reference  points  to  define  the  toroid.   Some  individuals 
may  take  exception  to  that  statement,  but  here  is  where  it 
comes  into  play.   Another  way  to  look  at  p,  the  radius  of 
curvature  in  the  4>   direction,  is  that  p  is  the  length  of  the 
arm  joining  the  point  in  question  Csay  point  P)  and  the  axis 
of  symmetry  Cthe  Z  axis  in  our  case)  and  coincident  with  the 
local  unit  normal  vector  at  point  P.   The  angle  6  can  be 
measured  either  at  the  axis  of  rotation  of  the  small  circle 
forming  the  shell,  or  at  the  intersection  of  the  line  p  and 
the  axis  of  symmetry  for  the  entire  toroid.   Now  when  6  ■  50° 
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for  instance,  the  geometry  looks  like  that  shown  in  Figure 
16.   As  e  decreases  toward  0,  tRe  end  of  p  must  travel 
farther  and  farther  down  the  axis  of  symmetry.   All  this  time 
p  is  getting  larger  and  x^   is  getting  smaller.   When  0   gets 
to  zero  plus  C0+)  the  situation  is  one  of  approximations.   p 
is  the  hypotenuse  of  the  right  triangle  with  the  right-angle 
at  the  origin.   Sin6  =  £  from  simple  trigonometry,  but  as  p 
goes  to  oo,  sinO  goes  to  zero  and  hence  6  must  go  to  zero 
also. 

Continuing  6  around  to  0",  the  closest  point  on  the  axis 
of  symmetry  in  line  with  the  normal  vector  is  now  at  +oo. 
This  now  gives  p  a  direction  opposite  to  that  of  the  local 
normal  vector,  hence  a  minus  sign,  and  as  6  goes  further 
negative,  p  is  now  smaller  in  magnitude.   CSimilar  arguments 
can  be  made  for  the  bottom  crown  where  6  *  180°)   The  sole 
reason  for  the  ±  <x>  trip  in  p,  and  hence  the  sign  change  in 
* .,    is  the  offset  R.   Somehow  this  must  be  related  to  the 
difficulty  in  solving  the  complete  toroid. 
GEOMETRY 

The  other  subtility  of  the  toroid  is  the  way  the 
mathematical  representation  of  the  surface  impacted  on  the 
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attempt  to  integrate  for  energy  over  the  toroid  surface. 

This  work  specifically  avoided  some  problems  encountered  by 

other  investigators,  significantly  the  expression  for  the 

curvature  in  the  <p   direction.   The  early  development  in  this 

work  specifically  avoided  the  term  C — - — ),  selecting  instead 

sin© 

an  expression  containing  the  term  Ca  +  sin©).   Force  balance 
and  constitutive  equations  were  developed  using  the  latter 
expression.   When  the  constitutive  equations  were  plugged 
into  the  energy  expression  CV),  no  singularities  existed. 
This  was  good  news. 

Unfortunately,  the  Ca  +  sin©)  term  also  appeared  in  the 
expression  lor  the  differential  surface  area  CdS).   The  bad 
news  was  that  many  of  the  elements  in  V  contained 
combinations  of  trigonometric  functions  not  readily  found  in 
the  CRC  Math  Tables  [27].   This  is  what  forced  the  numerical 
integration  of  the  elements  in  V. 

This  heavy  dependence  of  the  trigonometric  functions  in 
describing  the  problem  has  an  astounding  impact  on  the 
selection  of  an  assumed  functional  for  the  energy  method 
solution.   Any  combination  of  CA  ©  sin©)  or  C©A  sin©)  goes  to 
zero  when  integrating  from  0  to  2n.   Remember  that  the  dS 
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term  contains  a  sin©  term  so  all  integrals  contain  at  least 
one  sin©  term. 

It  should  be  pointed  out  that  the  solution  to  a 
symmetric  problem  can  be  obtained  by  solving  only  one  of  the 
unique  parts  of  the  problem.   This  was  not  done  in  this  case 
because  it  was  hoped  that  the  program  being  developed  would 
be  more  general  in  nature,  actually  leading  to  orthotropic 
analysis  of  a  stiffened  toroid. 

The  only  redeeming  quality  of  the  integrals  also  comes 
from  geometry.   Geometry  terms  containing  cos©,  cos2©,  sin©, 
sin2©,  and  (a  +  sin©>  and  thier  products  appear  in  numerators 
and  denominators  through  out  the  integrals.   Mister  Simpson 
and  the  digital  computer  were  allowed  to  seperate  the  wheat 
from  the  chaff. 

An  attempt  was  made  to  develope  a  non-trigonometric 
functional  which  would  satisfy  the  conditions  of 

connectivity  stated  earlier  and  to  meet  connectivity  at 
©  «  0°/360°.   The  conditions  to  be  met  were: 

f<0)  ■  fC360) 

f'<0)  =  f  <360> 

fC90)  ■  0 
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f  C90)  ■  0 

fC270>  -  0 

f'<270>  -  0 
This  resulted  in  eight  conditions.   A  functional  of  the 
following  form  was  generated: 

f<me>  =  A  +  BCmO>  +  CCme)2  +  DCm0)3 

+  ECmO)4  +  FCme>5  +  GCm6)6  +  HCmO)7 
A  was  arbitrarily  set  to  1  Cto  force  a  non  zero  solution  at 
6  =  O0).   Unfortunately,  by  the  time  fCmq)  was  processed, 
the  functional  looked  the  same  for  all  values  of  m.   The 
coefficients  all  changed  as  m  changed,  but  the  shape  of  the 
functional  remained  unchanged  from  one  value  of  m  to  the 
next  because  of  the  eight  conditions  imposed  upon  the 
functional.   This  lost  the  variety  of  the  solution  that  one 
expects  to  obtain  using  a  varying  m.   See  Figure  17. 

It  is  obvious  that  more  complex  functionals  are  needed. 
This  work  included  an  initial  evaluation  of  several 
functions  for  this  task.   A  considerable  amount  of  time  was 
put  into  investigating  the  hypergeometric  functions, 
particularly  the  Legendre  functions.   In  addition  to  serving 
as  a  functional  for  u  and  w,  it  was  also  hoped  that  this 
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function  would  satisfy  the  form  of  the  differential  equation 
for  the  toroid  and  hence  yield  an  exact  solution.   (It  was 
partially  the  anticipation  of  using  such  a  complex  function 
that  drove  the  layout  of  the  program  Csee  Appendix)  to  have  a 
seperate  generating  section  for  the  functional.)   The 
Legendre  function,  or  any  other  hypergeometric  function,  was 
not  used  due  to  compatability  conditions  being  harder  to 
achieve.   The  next  step  in  this  progression  would  be  the 
assumption  of  a  combined  trigonometric  function  for  u  and  w, 
Cie.  A  sinCmO)  +  B   cosCmO)). 
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CLOSING 

In  Chapter  1  it  was  stated  that  the  geometric  curvature 

and  the  impact  of  shell  geometry  on  the  assumed  functional 
were  the  two  problem  areas  in  this  thesis.   Of  the  two,  the 
former  is  the  more  significant  in  this  writers  opinion.   The 
analysis  behind  Figure  15  is  crude,  but  it  is  an  attempt  to 
represent  in  graphics  and  mathematics  what  the  shell  is 
doing.   Being  unable  to  do  this  is  what  makes  the  solution  so 
difficult. 

In  a  dynamics  class,  one  student  asked  the  professor  if 
another  mathematical  technique  could  be  used  to  analyse  a 
spring.  The  professor  thought  briefly,  then  answered,  "You 
can  do  all  the  math  you  want  but  you  have  to  think  like  a 
spring.  The  spring  knows  what  it  is  doing. "  This  writer 
feels  a  little  short  on  thinking  like  a  toroid,  but  also 
feels  in  good  company. 
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APPENDIX 


THE  PROGRAM 


This  program  is  to  set  up  the  geometry,  the  energy 
expression.  Simpson's  multipliers,  and  perform  the  matrix 
manipulations  to  solve  for  the  displacements  <u>  and  <w>. 


50  CLS: PRINT: PRINT: INPUT  "WHAT  OUTPUT  FILE  NAME  "jAS 
60  CLS 
1000  NN=  50 


1010  MM=10 


1020 
1030 
1040 
1050 
1060 
1070 
1080 
1090 

1100 
1110 
1120 
1130 
1140 
1150 
1160 
1170 
1180 
1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1290 
1300 
1310 
1320 
1330 
1340 
1350 
1360 

1370 
1380 
1390 
1400 
1410 
1420 
1430 
1450 


rsu=.  j 

R=2 

A=4 

H=.02 

E=3E+07 

P=-100 


":A*R;"  inches" 
":R: "  inches" 
";H;"  inches" 


:NU 
:E 


*This  is  the  number  of  steps  to  be  used  in  the 
Simpsons  integral 

*This  is  the  number  of  terms  in  the  series  for 
displacements  u  ans  w. 

*Poisson"s  ratio 

*R  from  the  text. 

*a  from  the  text. 

*t  from  the  text. 

♦Young's  modulus 

♦pressure  in  psi. 
CLS: PRINT: PRINT:  PRINT 
PRINT  "The  following  are  the  default  values  in  the 

program: " 
PRINT: PRINT 
PRINT  "GEOMETRY" 
PRINT  "Radius  of  rotation 
PRINT  "Radius  of  the  circle 
PRINT  "Thickness  of  the  shell 
PRINT: PRINT 
PRINT  "MATERIAL" 
PRINT  "Poisson's  ratio 
PRINT  "Modulous  of  elasticity 
PRINT: PRINT 
PRINT  "SOLUTION- 
PRINT  "Number  of  steps  in  the  Simpsons" 
PRINT  "   integral 

PRINT  "Number  of  terms  to  be  used  in" 
PRINT  "   the  series  approximation 
PRINT  "Pressure 
PRINT: PRINT 

INPUT  "Do  you  wish  to  change  any  of  these  CY/TO 
IF  SS="N"  THEN  1520 
CLS 

PRINT: PRINT 
PRINT  "GEOMETRY- 
INPUT  "Radius  of  rotation  (inches) 
INPUT  "Radius  of  the  circle  Cinches) 
A=RR/R 

IF  A  >  1  THEN  1380 
CLS:PRINT"RADIUS  OF  ROTATION  MUST  BE  GREATER  THAN  THE 

RADIUS- 
PRINT  "'OF  THE  CIRCLE.  ":  GOTO  1300 
INPUT  "Thickness  of  the  shell  Cinches) 
PRINT: PRINT 
PRINT  "MATERIAL" 
INPUT  "Poisson  s  ratio 
INPUT  "Modulus  of  elasticity 
PRINT: PRINT1440  PRINT  "SOLUTION" 
PRINT  "Number  of  steps  in  the  Simpsons" 


";NN 

";MM 
";P;"  psi" 


":S$ 


";RR 


";H 


";NU 
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1460  INPUT  "   integral  : M;NN 

1470  PRINT  "Number  of  terms  to  be  used  in*' 
1480  INPUT  "   the  series  approximation     :":MM 
1490  INPUT  "Pressure  (psi)  :";P 

1500  ' 
1510  ' 

1520  'THIS  SECTION  IS  TO  SET  UP  THE  PARAMETERS  FOR  THE 

INTEGRATION  SECTION 
1530  ' 

1540  DT=6.283185307#/NN 
1550  I=NN+1 

1560  DIM  PC  I.  3,  MM)   ^Defines  the  functions  <p. 
1570  DIM  JCI)       *Simpson's  multiplier 

Arrays  MD,  ME,  MF,  and  MG  correspond  to  the  matrixes  in 
the  partial  of  V  with  respect  to  u,  and  w,.   MP  is  the  matrix 

representing  the  energy  due  to  the  pressure  term.   MS  is  a 
scratch  matrix.  M*^  is  a  series  of  matrixes  used  in  the 
matrix  manipulation  subroutines.   MV  and  MU  are  the  answers. 

1580  DIM  MD(MM.MM):DIM  MECMM, MM) : DIM  MF(MM,MM) 

1590  DIM  MG(MM.MM):DIM  MS(MM,MM) 

1600  DIM  MX*(MM.MN):DIM  MY#CMM, MM) : DIM  MZ#CMM,MM) 

:DIM  MI^CMM.MM) 
1610  DIM  MP(MM):DIM  MV(MM) : DIM  MU(MM) 
1620  CLS: PRINT: PRINT: PRINT  "GENERATING  FUNCTIONS" 
1630  FOR  I  =  1  TO  NN+1 
1640   T=CI-1)*DT 
1650   FOR  M  =  1  TO  MM 

1660   PCI.1,M)=SIN(M*T)  'This  is  the  function 

1670   P(I.2,M)=M*C0S(M*T)        'This  is  the  first 

derivative 
168-0   P(I.3,M)=-M'2*SIN(M*T)     'This  is  the  second 

derivative 
1690   NEXT  M 
1700   J(I)=4 

1710   IF  (l/2-INT(I/2))>0  THEN  J(I)=2 
1720  NEXT  I 
1730  J(l)=l: J(NN+1)=1 

Integrate  using  Simpson's  Rule 

2000  'INTEGRATE  FROM  0  TO  2*PI 

2010  IC=3. 14159*(E*H/(l-NU~2))*DT/3 

2020  BC=H'2/(6*R'3) 

2030  PC=-2*3.14159*P*R'2*DT/3 

2040  FOR  I  =  1  TO  NN+1 

2050   T=(I-1)*DT 

2060   S=SIN(T) 

2070   C=C0S(T) 

2080   CLS: PRINT: PRINT: PRINT"INTEGRATING  ON  STEP  ":I 

2090   PRINT  "  ANGLE  ":T*57.296;"  DEGREES" 

2100   AA=A+S2102   HH=H' 2/(12*R~2) 

2104   CC=C/AA 

2106   CA=2*NU*CC 

2108   CB=CC2 

2110   CD=2*NU*S/AA 

2112   AX=2*S*C/AA'  2 
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2200 


2130  FOR  K  =  1  TO  MM 

2140  MP(k)=MP(k)-J(i;>*P(I.l,k;>*AA 

2150  SQ=1/AA'2 

2160  TR=2*NU/AA 

2170  B\=H'3*Q+2*NU*C/AA+CC/AA)~2)/a2*ir2) 

2180  FOR  M   =  1    TO  MM 

2190        MDCk.M)=MD(k,M)+ja)*(a+HH)*(2*PCI,3,K)*P(I,3,M) 

+CA*P(I.3,1C)*PCI,2,M$ 
+CB*2*pa,2,K)»HP(Ii2rM)))*AA 
ME(K,M)=ME(K,M)+J(I)*(2*f>(t ,  i^M)*£cL  3, 1C) 

+CA*(P(  1 .  3^*?(  1 . 1  Ji$*& 
+P(I,liM)*f>(I^2iK>*6) 
+AX*p(L2<K)*£cijL<M) 
-HH*(2*f>(L  3^M)*f>(t   3^K) 

+CA*(f>ct,3^K)*f>(L2,M) 
+PCI.3.m5*PCI.2.k5) 
+2*CB*P( I , 2 . K) *P( I . 2 , M) ) ) *AA 
2210        MF(k.M)=MFCk.M)+J(I)*(2*PCI,3,M:>*PCI.l,K5 

+2*NU*CP(  1 . 1 .  K)*P(  1.3.  M)*S 
+P(Ii2iM)*PClil^K)*C)/AA 
+AX*P(I^2^M)*P(t.l.K) 
-HH*C2*P( I . 3 . M)*P( I . 3^K) 

+CA*(P(L32M)*P(t;2,D 
+P(I.3iK3*PCI,2iM$) 
+2*CB*£(I,2,M}*kl,2,K;:>:>:>*AA 
2220        MGCk.M)=MG(k,M)+J(I)*(2*P<:i,l,M5*P<:i,l,K.5 

+2*CD*P(  I   1  r  f|)*PC  I,  1 ,  fo 
+2*P( I . 1 . M)*P( I . 1 , K)*S~2/AA~2 
+HH*C2*P<:iz3.M)*P(I.3iK) 

+CA*CP(I.3iM)»Hp(i.2,K.) 

+Pa,3JC>*Pa.2.M5) 

+2*CB*PCI,2,M5^CI,2,K)))*AA 


2230 

NEXT  M 

2240 

NEXT  K 

2250 

NEXT  I 

2260 

> 

2270 

y 

2280 

FOR  K  =  1  TO  MM 

2290 

MP(k)=POMPCk) 

2300 

FOR  M  =  1  TO  MM 

2310 

MD(k,M)=IC*MD(k, 

M) 

2320 

ME(k,M)=IC*ME(k, 

M) 

2330 

MF(k.M)=IC*MF(K. 

M) 

2340 

MG(k.M)  =  IOMG(k. 

M) 

2350 

NEXT  M 

2360 

NEXT  k 

3000  'SOLVE  THE  MATRICES  FOR  w  AND  u 

m        m 

3010  CLS: PRINT: PRINT: PRINT  "SOLVING" 

3020  ' 

3030  FOR  k  =  1  TO  MM:  FOR  M  =  1  TO  MM: MX#CR, M)=MD(k, M) : NEXT 

M:NEXT  k 

3040  GOSUB  10000  '  He****  Invert  MD  **** 

3050  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM 

3060   MY^(K,M)=MECR.M):MX#(K,M)=-MI#(K,M) 

3070  NEXT  M:NEXT  R 

3080  GOSUB  11000  '  ****  (-MD'-1)*ME  **** 

3090  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM 

3100   MY*(k.M)=MZ^Ck.M):MX#(k.M)=MFCk,M) 

3110  MS(k.M)=MZ*(k.M)  '      ****  Saving  this  for  Urn 

solution  **** 

3120  NEXT  M:NEXT  K. 
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3130  GOSUB  11000  '  ****  MF*C-MD' -1)*ME  **** 

3140  FOR  K  =  1  TO  MM: FOR  M  =  1  TO  MM 

3150   MY#CK, M)=MGCk. M)  :  MX#Ck, M)=MZ#CK, M) 

3160  NEXT  M:NEXT  H 

3170  GOSUB  12000  '  ****  MG  -  MF*C-MD~-1)*ME 

3180  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM: MX#CK. M)=MZ#(K, M) 

:NEXT  M:NEXT  K 
3190  ' 

3200  'ROV  REDUCTION 
3210  'Input  is  MX#CK.M)  and  MPCK) 
3220  'Output  is  the  altered  MX^CK.M)  and  MPCK) 
3230  FOR  S  =  1  TO  MM-1 
3240   FOR  K  =  S+l  TO  MM 
3250   Q=-MX#(K.S)/MX^(S.S) 
3260   MPCK)=Q*MPCS)+MPCK) 
3270   FOR  M  =  1  TO  MM 
3280    MX,*CK.M)=Q*MX#CS.M)+MX#CK.M) 
3290   NEXT  M 
3300   NEXT  K 
3310  NEXT  S 
3320  ' 

3330  'SOLVE  FOR  MVCM) 

3340  'Continuing  on  using  MX#Ck.M)  and  MPCK)  from  before 

3350  FOR  S  =  MM  TO  1  STEP  -1 

3360   MV/CS)=MPCS) 

3370   FOR  M  =  S+l  TO  MM 

3380   MVCS)=MVCS)-MVCM)*MX,*CS.M) 

3390   NEXT  M 

3400   MWCS)=MVCS)/MX#( S . S) 

3410  NEXT  S 

3420  " 

3430  '  SOLVE  FOR  MUCK) 

3440  FOR  K  =  1  TO  MM 

3450   FOR  M  =  1  TO  MM 

3460    MUCk)=MUO;)+MSCk.M)*MVCM) 

3470   NEXT  M 

3480  NEXT  k 

3490  '  PRINT  OUT  RESULTS 

3500  PRINT  T,"  Vm  ","  Urn  " 

3510  FOR  M  =  1  TO  MM 

3520   PRINT  M. MVCM). MUCM): PRINT 

3530  NEXT  M 

3540  OPEN  "0".#1,AS   *This  section  puts  the  output  onto 

3550  PRINT  *1."NN  =  ";NN:"  '    MM  =  ";MM 

3555  PRINT  #1 '.  "RADIUS  =";A*R;"      radius  =  ";R 

3560  PRINT  #1,"  " 

3570  PRINT  #1.  "#","   vm  "."  Urn  " 

3580  FOR  M  =  1  TO  MM 

3590   PRINT  n,    M,MVCM),MUCM):PRINT 

3600  NEXT  M 

3610  CLOSE 

3620  END 
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This  section  contains  routines  used  in  the  main  program. 

10000  'INVERSION  OF  A  MATRIX 

10010  'Input  is  MX#(K,M) 

10020  'Output  is  MI#(LM) 

10030  FOR  K  =  1  TO  MM 

10040   MI#CK.K)=1 

10050  NEXT  K 

10060  FOR  S  =  1  TO  MM-1 

10070   FOR  K  =  S+l  TO  MM 

10080   Q=-MX#(K,S)/MX#CS.S) 

10090   FOR  M  =  i  TO  MM 

10100    MX#(K,M)=0*MX^(S.M)+MX#(k.M) 

lono   mi#ck.m)=q*mi*cs.m)+mi#(k.m) 

10120   NEXT  M 

10130   NEXT  K 

10140  NEXT  S 

10150  FOR  S  =  MM  TO  1  STEP  -1 

10160   FOR  R  =  S-i  TO  1  STEP  -1 

10170   Q=-MX#(K.S)/MX^CS^S) 

10180   FOR  M=MM  TO  1  STEP  -1 

10190    MX#CK.M)=Q*MX#CS,M)+MX#CK.M> 

10200    MI*CK'.M>=0*MI#CS.M)+MI#CK'.M> 

10210   NEXT  M 

10220   NEXT  K 

10230  NEXT  S 

10240  FOR  K  =  1  TO  MM 

10250   FOR  M  =  1  TO  MM 

10260   MI#CK,M)=MI#CK.,M>/MX#CK,K) 

10270   NEXT  M 

10280  NEXT  K 

10290  RETURN 

11000  'MULTIPLY  TWO  MATRICES 

11010  'Multiply  MX^CK.M)  times  MY#(K,M)  in  that  order! 

11020  'The  output  is  MZ^Ck.M) 

11030  FOR  K  =  1  TO  MM 

11040  FOR  M  =  1  TO  MM11050   MZ#(K.M)=0 

11060  FOR  S  =  1  TO  MM 

11070  MZ^(K.M)=MZ^(R,M)+MX^(R.S)*MY#(S^M) 

11075  IF  ABSCMZ#CK,M$X.  0000005  THEN  MZ#(K.M)=0 

11080  NEXT  S 

11090  NEXT  M 

11100  NEXT  K 

11110  RETURN 

12000  'ADD  TV/0  MATRICIES 

12010  'Add  MX^CK.M)  to  MY^(K.M) 

12020  'Output  is  MZ#CK,M) 

12030  FOR  K  =  1  TO  MM 

12040   FOR  M  =  1  TO  MM 

12050   MZ^(K.M)=MX#CR,M)+MY#(K,M) 

12060   NEXT  M 

12070  NEXT  K. 

12080  RETURN 
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